%
%一个一般系统的哈密顿量
%
function [grn, grn_t1, grn_t2] = grn_gnrl(hmlt, temp, nt, dtau)
    hsize = size(hmlt);
    %grn = zeros(hsize(1), hsize(2));
    [vecs, vals] = eig(hmlt);
    weight_mat = zeros(hsize(1), hsize(2));
    if temp < 0.01
        for idx = 1:hsize(1)
            if vals(idx, idx) < 0
                weight_mat(idx, idx) = 1.0;
            elseif vals(idx, idx) == 0
                weight_mat(idx, idx) = 0.5;
            else
                weight_mat(idx, idx) = 0.0;
            end
        end
        allow = sum(sum(weight_mat))
    else
        for idx = 1:hsize(1)
            weight_mat(idx, idx) = 1.0 / (1.0 + exp(vals(idx, idx)/temp));
        end
    end
    grn = vecs * weight_mat * conj(transpose(vecs));
    if temp < 0.01
        warning('only equal time green function is implement')
        grn_t1 = 0;
        grn_t2 = 0;
        return
    end
    %grn_t1 = c_i(tau) c^d_j(0) = 
    %e^{tau H} c_i e^{-tau H} c^d_{j} =
    %sum_{k}e^{tau H} <k|i> c_k e^{-tau H}  <j|k> c^d_{k}
    % sum_{k} e^{-tau E_k} <k|i>  1/(1+e^{-beta E_k})  <j|k>
    grn_t1 = zeros(hsize(1), hsize(2), nt+1);
    for tauidx=0:1:nt
        tau = tauidx * dtau;
        wmat_t1 = zeros(hsize(1), hsize(2));
        for idx = 1:hsize(1)
            wmat_t1(idx, idx) = exp(-tau*vals(idx, idx)) / (1.0 + exp(-vals(idx, idx)/temp));
        end
        grn_t1(:, :, tauidx+1) = vecs * wmat_t1 * conj(transpose(vecs));
    end
    %grn_t2 = c^d_i(tau) c_j(0) = 
    %e^{tau H} c^d_i e^{-tau H} c_{j} =
    %sum_{k}e^{tau H} <i|k> c^d_k e^{-tau H}  <k|j> c_k
    % sum_{k} e^{tau E_k} <i|k>  1/(1+e^{beta E_k})  <k|j>
    grn_t2 = zeros(hsize(1), hsize(2), nt+1);
    for tauidx=0:1:nt
        tau = tauidx * dtau;
        wmat_t2 = zeros(hsize(1), hsize(2));
        for idx = 1:hsize(1)
            wmat_t2(idx, idx) = exp(tau*vals(idx, idx)) / (1.0 + exp(vals(idx, idx)/temp));
        end
        grn_t2(:, :, tauidx+1) = vecs * wmat_t2 * conj(transpose(vecs));
    end
end

